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Employing smart Monte Carlo sampling techniques within the grand canonical ensemble, we 
investigate the properties of water at a model hydrophobic substrate. By reducing the strength 
of substrate-water attraction we find that fluctuations in the local number density, quantified by 
a rigorous definition of the local compressibility x( z )i increase rapidly for distances z within 1 
or 2 molecular diameters from the substrate as the degree of hydrophobicity, measured by the 
macroscopic contact angle 9, increases. Our simulations provide evidence for a continuous (critical) 
drying transition as the substrate-water interaction becomes very weak: cos (9) —> —1. We speculate 
that the existence of such a transition might account for earlier simulation observations of strongly 
enhanced density fluctuations. 


All physical scientists would agree that for water at a 
flat substrate a contact angle 6 > 90° defines the sub¬ 
strate as hydrophobic. For such values Young’s equa¬ 
tion implies that the substrate-vapor interfacial tension is 
lower than the substrate-liquid tension, i.e. the substrate 
prefers vapor to liquid. 9 is determined by surface chem¬ 
istry, which in turn determines the strength and range of 
substrate-fluid interactions. A key question, much dis¬ 
cussed by the chemical physics communities, is whether 
there is an effective indicator of local ordering of wa¬ 
ter, manifest at microscopic distances from the substrate, 
which might correlate with the degree of hydrophobicity 
as measured by the macroscopic thermodynamic quantity 
9. Quantifying the character and spatial extent of water 
ordering at hydrophobic entities is important across sev¬ 
eral disciplines, ranging from applied physics and mate¬ 
rials science, where understanding slip lengths of water, 
shown to be correlated with substrate hydrophobicity, is 
crucial in microfluidics [1], to bio-physical processes such 
as protein-folding and micelle and membrane formation 
[2]. Many studies have focussed on the average one-body 
density, i.e. the density profile of oxygen atoms near 
the substrate. Experimentally this is difficult to measure 
and several conflicting results were reported. Neverthe¬ 
less, by 2009 a consensus emerged from x-ray reflectiv¬ 
ity measurements that for water near a variety of hy¬ 
drophobic self-assembled-monolayers (SAMs), there is a 
region of density depletion corresponding to only a frac¬ 
tion of a water monolayer [3-5]. This was challenged by 
Chattopadhyay et.al. [6], who reported larger depletion 
lengths, increasing with 9, for water at fluoroalkylsilane 
SAMs. For the largest contact angle 9 = 120° the deple¬ 
tion length was 8A, corresponding to about three water 
layers. However, in a re-analysis [7] of the x-ray reflec¬ 
tivity data of Ref. [6] it was argued that such a large 
hydrophobic gap is likely to be an artefact of the data 
analysis that was used. This is disputed [8]. Certainly 
such large depletion lengths are at odds with results of 
many computer simulations of water models, where the 


thickness of the depleted density region varies typically 
between 1.5 — 2.0A for 9 between about 110° and 130° 
[9], 

Other simulation studies have focussed on the density 
fluctuations of water at model hydrophobic substrates, 
arguing that some measure of the local compressibility 
might provide a better indicator of hydrophobicity than 
does the local density profile. The group of Garde [10] 
takes this stance and the review by Janradagni et. al. 
[11] places their work in the context of bio-molecular sys¬ 
tems. Chandler and co-workers [12-14] and Mittal and 
Hummer [15] also emphasize that density fluctuations 
can be enhanced at hydrophobic substrates. Although 
these studies identify some underlying phenomenology, 
the various measures of the local compressibility intro¬ 
duced are ad-hoc. In particular these do not correspond 
to integrals over density correlation functions in the in¬ 
homogeneous liquid. Clear insights into the nature of 
the underlying fluctuations requires such a measure [16]. 
Choice of ensemble is important. Most members of the 
community simulating water, choose not to work grand 
canonically. In most real interfacial situations there is 
a reservoir and it is natural to vary the chemical poten¬ 
tial p of the liquid. Following a recent analysis [17] of 
density fluctuations in a Lennard-Jones (LJ) liquid at a 
planar substrate (wall), we argue that the most appropri¬ 
ate definition [16] of the local compressibility, for a fixed 
confining volume, is 

X(z) = (>dp(z)/dp) T , (1) 

where p(z) is the average one-body density, z is the 
distance normal to the substrate, and the temperature 
T is fixed. For a bulk fluid with constant density pi, 
x(z) = Xb = Pb K T, where kt is the isothermal compress¬ 
ibility. The definition (1) is consistent with x( z ) as an 
integral over the density-density correlation function and 
yields a unique fluctuation formula for this quantity; See 
Eqs. (8,9) of [16]. In ref. [17] it was shown using classi¬ 
cal density functional theory (DFT) that for 2 : close to 



2 


the substrate x( z )/Xb becomes large as the substrate be¬ 
comes more solvophobic, i.e. 9 increases. The effects are 
not small: for 9 ~ 160° the ratio is about 25. 

In this Letter we determine the local compressibility of 
the extended simple point charge (SPC/E) model of wa¬ 
ter [18] near a planar substrate. On reducing the strength 
of substrate-water attraction, thereby increasing the hy- 
clrophobicity as measured by 9, we find x( z ) increases in 
a similar fashion to the LJ case [17]. We focus on the ap¬ 
proach to complete drying, cos (9) —> — 1 at bulk vapor- 
liquid coexistence /i = p cx . Determining the nature of 
this transition continues to be challenging, even for sim¬ 
ple fluids at solvophobic substrates. Using smart sam¬ 
pling techniques within Grand Canonical Monte Carlo 
(GCMC) simulations we find evidence for a critical dry¬ 
ing transition in SPC/E water at a weakly attractive sub¬ 
strate. We note that recent MD simulations [14, 19] of the 
same water model also investigate very weak substrate- 
water attraction. Although both studies attempt to link 
the growth of density fluctuations to increasing 9, nei¬ 
ther simulation measures 9 accurately so proximity to 
the drying point is uncertain. Moreover, neither makes 
explicit the possibility of a critical drying transition. We 
argue the latter is key to understanding the properties 
of water models in the extreme hydrophobic regime and, 
in particular, the enhanced density fluctuations that are 
observed. 

We choose SPC/E as a simple but realistic model of 
water, known to provide a reasonable account of bulk 
vapor-liquid coexistence and the vapor-liquid surface ten¬ 
sion of real water [20, 21]. Moreover, this model has of¬ 
ten been employed in simulations of water at hydropho¬ 
bic substrates; for examples of MD studies see refs. 
[10, 14, 19, 22] and for GCMC studies see refs. [20, 23- 
25]. Kumar and Errington employ simulation techniques 
similar to the present but they do not measure xi z )> i- e - 
they do not access local density fluctuations. Rather they 
measure [23] the surface excess compressibility Xex which 
is an integrated measure of ‘excess’ density fluctuations 
throughout the system [16, 17, 26] 

Our simulations were performed at T = 298K. Since 
liquid water at room temperature is too dense for stan¬ 
dard GCMC to operate effectively, smart sampling tech¬ 
niques were implemented. Configurational Bias MC was 
used to insert, delete, translate and rotate molecules [27], 
while Transition Matrix MC [28], Multicanonical Sam¬ 
pling [29] and Histogram Reweighting [30] were deployed 
to smoothly connect the vapor and liquid regions of con¬ 
figuration space. Together these methods allowed us to 
simulate modest system sizes of order a few hundred wa¬ 
ter molecules very accurately. Although MD simulations 
typically deal with greater numbers of molecules, the 
GCE is the appropriate ensemble for accurately studying 
x(z) and thermodynamic quantities such as cos(f?) and 
surface phase behaviour because it permits direct study 
of the fluctuations which characterise phase transitions. 


As these fluctuations occur on the scale of the system 
itself, the GCE is less afflicted by finite-size effects than 
other simulation ensembles. 

We consider two simulation setups: (i) a fully peri¬ 
odic cubic box of volume V = L 3 ; (ii) a semi-periodic 
cuboidal slit geometry of volume V = L 2 D, in which 
the oxygen atoms interact with a pair of symmetry- 
breaking walls separated by distance D. For the latter 
geometry, the single wall-oxygen potential takes the form 
V s (z) = ae w f[2/15{a w f/z) 9 -(a w f/z) 3 ], where e wf is the 
wall-fluid interaction strength and 2 is the distance from 
the wall, a is a constant, the choice which sets the units 
of energy (see Supplementary Material [16]), while a w f 
sets the length scale for wall-fluid interactions which we 
assigned to be 3.5A, in accordance with previous studies 
[23]. 

The choice of e w f controls the degree of hydrophobic- 
ity, and hence the contact angle 9. Our use of the GCE 
permits us to calculate 9 directly from Young’s equation, 

Xvi cos (9) — ^fwv Xwi • (2) 

Here 7 v i is the vapor-liquid interfacial tension and j wv 
and 7 w i are the wall-vapor and wall-liquid interfacial ten¬ 
sions, respectively. Provided one can calculate "f v i and 
lwv~lwi at n = n cx , cos (9) is given. Both quantities are 
directly obtainable from measurements of the probabil¬ 
ity distribution p{p) of the fluctuating molecular number 
density p = N/V in the appropriate simulation geometry. 
Specifically, studies of p{p) in the fully periodic system 
allow estimates of both p cx and 7 v i , while studies of p(p) 
in the slit geometry allow measurements of j wv — 7 w i. 

For the periodic system at vapor-liquid coexistence, 
p(p) exhibits a pair of equally weighted peaks [31] sep¬ 
arated by a flat probability ‘valley’. The low and high 
density peaks correspond to pure vapor and liquid phase 
states respectively, while the flat valley corresponds to 
liquid-slab configurations in which a pair of liquid-vapor 
interfaces align parallel to one face of the simulation box. 
Since in the probability valley p m in is typically many 
decades smaller than the peak probabilities p v ap and pn q , 
standard sampling cannot plumb the valley depth. For 
this reason, Transition Matrix and Multicanoncial MC 
techniques [28] were used to accumulate p{p) in histogram 
form across the full range of density from vapor to liquid. 
Initially p(p) was determined for a near-coexistence state 
point. The distribution was then reweighted [30] with 
respect to p to determine p cx via the equal peak weight 
criterion. Once obtained, the coexistence form of p(p) 
permits an estimate of 7 v i [32, 33]: 

7 vi = (2/3 L 2 )- 1 In(p max /p min ) , (3) 

where /3 = (fc B T ) _1 and p max = p vap = Piiq is the height 
of the pure phase peaks. 

In a similar manner, estimates of p{p) at coexistence 
were accumulated for the slit geometry at various e w f. 
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The resulting histograms (Fig. 1) typically exhibit two 
peaks -one at low density corresponding to wall-vapor 
(i.e. vapor at the wall) configurations, and another at 
high density corresponding to wall-liquid (liquid at the 
wall) configurations [34], For a given e w f, 7 wv — "f w i is 
calculated from the measured ratio of vapor and liquid 
peak heights in p(p) [35]: 

7 wv 7 wi = (2/3T“) ln(p va p/piiq) . (4) 



FIG. 1. (Color online), (a) The measured forms of p(p) for 
a selection of values of e w f spanning the range from wetting 
to drying; the contact angles 9 are shown. The system size is 
L = D = 20A. (b) A closeup of the region close to the drying 
transition. Note the continuous erosion of the liquid peak as 
e w f decreases. In (a) and (b) the dashed line is p(p) for the 
periodic system with L = 20A. 

We now discuss the pertinent features of Fig. 1, not¬ 
ing that —In p{p) measures the grand potential of the 
fluid. For large values of e w f, wall-liquid configurations 
are much more probable than wall-vapor configurations. 
Indeed for e w f ~ 10.3, we find cos(0) = 1, correspond¬ 
ing to the transition to complete wetting by liquid wa¬ 
ter. Reducing e u ,f below this value takes the system 
first into the partially wet (hydrophilic) regime for which 
0° < 0 < 90°. Thereafter, for e w f < 6 , the system 
enters the partially dry (hydrophobic) regime in which 


90° < 6 < 180° and where wall-vapor configurations are 
favoured over wall-liquid ones [36]. 

As e w f is reduced within the hydrophobic regime, 
e w f < 6 , Fig. 1 shows that the height of the liquid peak 
diminishes progressively, until it vanishes smoothly into 
a plateau in In p(p) at a low wall strength which, for this 
system size, is e w f = 0.5(1). Interestingly, the vanishing 
of the liquid peak occurs precisely at the drying point 
cos(0) = —1. This is clear from the dashed curve in 
Fig. 1(b) which shows the form of In p(p) corresponding 
to the fully periodic system at coexistence, for the same 
value of L. The peak to valley separation in this plot is 
the same as the vapor to liquid peak separation in the slit 
system at e w f = 0.5. It follows from Eqs.(2)-(4) that the 
vanishing of the liquid peak in lnp(p) marks the drying 
point cos( 0 ) = — 1 . 

The behaviour of lnp(p) in the slit system as a func¬ 
tion of e w f reveals interesting qualitative differences be¬ 
tween the nature of the approach to complete wetting 
(cos( 0 ) —► 1 ), and to complete drying (cos (0) —► — 1 ). 
In the former case, the wall-vapor configurations remain 
metastable at the transition, as evidenced by the presence 
of a vapor peak at the wetting wall strength e w f = 10.3. 
This signifies that the wetting transition of water is a first 
order surface phase transition. By contrast, on approach¬ 
ing the drying point, the liquid peak vanishes smoothly 
into a plateau of constant probability. The distinction 
is borne out by a plot of cos (6) + 1 versus e w f calcu¬ 
lated from (2) and shown in Fig. 2. We find that cos(0) 
approaches unity with a non-zero gradient denoting un¬ 
ambiguously a first order wetting transition, but appears 
to approach —1 tangentially. Such a scenario implies a 
continuous (critical) transition to drying. 
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FIG. 2. Estimates of cos($) + 1 as a function of e w f spanning 
the region from wetting cos($) = 1 to drying cos($) = — 1 . 
Statistical uncertainties are smaller than symbol sizes. 

The physical interpretation of the smooth erosion of 
the liquid peak on the approach to drying is that the free 
energy barrier attaching the liquid to the wall vanishes 
continuously. This allows the liquid layer to detach and 
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be replaced by vapor. In our slit system the emergent 
liquid-vapor interface fluctuates freely (cf., the plateau in 
p(p) at e w f = 0.5) until the liquid slab thickness decreases 
sufficiently that the system undergoes (capillary) evapo¬ 
ration. This interpretation is confirmed by measurements 
of the number density profile p(z) of the oxygen atoms in 
the hydrophobic regime (Fig. 3). For medium strength 
attractive wall-fluid interactions which give rise to ‘neu¬ 
tral’ hydrophobicity, i.e. 6 ~ 90°, the liquid density is 
high at the wall and packing effects occur. The structure 
in the profiles decreases smoothly as e w f is reduced and 
a low density vapor layer starts to appear near the wall, 
the thickness of which appears to grow continuously as 
drying is approached. The thickness of the vapor layer 
represents the order parameter for the drying transition 

[37]. 



FIG. 3. (Color online). Normalized number density profiles 
p(z) for various values of e w f in the hydrophobic regime, pi, 
is the bulk liquid density at p cx - 

The behaviour of p(z) on the approach to drying re¬ 
veals the emergence of a vapor layer or density gap as¬ 
sociated with hydrophobic surfaces. However, a much 
more sensitive and revealing measure of the degree of hy¬ 
drophobicity is the compressibility profile xi z ) defined 
in (1), which provides a robust measure of local den¬ 
sity fluctuations close to the wall. Accurate estimates of 
this quantity are readily obtained within the GCE, by 
exploiting histogram reweighting to numerically differen¬ 
tiate the density profile p(z \p). The measured forms of 
x(z), normalized with respect to the bulk compressibil¬ 
ity Xb (i- e - the compressibility far from the wall), are 
presented in Fig. 4. These show that close to the wall, 
X{z)/Xb grows rapidly as e w f is reduced towards the dry¬ 
ing point, eventually exceeding its bulk value by nearly 
two orders of magnitude [38]. This finding mirrors what 
is found in DFT calculations of the solvophobic regime 
in a Lennard-Jones system [17], and reflects the develop¬ 
ment of a large transverse correlation length for density 
fluctuations as cos(0) —> —1 [16]. We note that there 


is a large enhancement in the local compressibility even 
for values of e w f which correspond to contact angles not 
much greater than 0 = 90° [39]. This confirms the rel¬ 
evance of our findings for experimental studies on real 
smooth hydrophobic surfaces such as SAMs where con¬ 
tact angles have a maximum value of about 130°. 



z (A) 


FIG. 4. (Color online). The normalized local compressibility 
X( x )/Xb for a selection of values of e w f. These correspond to 
the density profiles of Fig. (3). 

In summary, we have shown that smart sampling 
within the GCE is a powerful approach for accurately 
characterising the hydrophobic regime including the ap¬ 
proach to drying (6 —> 180°) in water. The local com¬ 
pressibility profile x( z ) is the proper statistical mechanics 
measure of local density fluctuations in fluids [16]. When 
applied to water in contact with a hydrophobic substrate, 
x{z) provides a sensitive indicator of how the micro¬ 
scopic structure near the substrate reflects the macro¬ 
scopic contact angle 6 -much more so than the density 
profile p(z) alone. In the drying limit the order param¬ 
eter (i.e. the thickness of the vapor layer) grows con¬ 
tinuously but slowly with decreasing e w f. In contrast 
x( z ) grows rapidly and exceeds its bulk value by nearly 
two orders of magnitude over the range of e w f explored, 
indicating the growth of a large transverse correlation 
length [16]. These findings point to drying in water at 
models of hydrophobic substrates as being a surface crit¬ 
ical phenomenon. Indeed separate simulation studies of 
drying in the Lennard-Jones system reproduce the phe¬ 
nomenology seen in our studies of water whilst increased 
system size permits the extraction of critical power law 
behaviour [40]. With these insights, one can rational¬ 
ize and explain the observations of enhanced fluctuations 
in previous simulation studies of water near hydropho¬ 
bic surfaces [10-15, 23]. We believe that these are at¬ 
tributable to the proximity of a surface critical point i.e. 
the approach to a continuous drying transition, the ef¬ 
fects of which extend throughout the hydrophobic regime 
but were not recognised previously. 

Some of the simulations described here were performed 
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on the Bath HPC Cluster. We thank R. Jack for helpful 
discussions. 
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